// 版权归Go作者所有。版权所有。
// 此源代码的使用受BSD样式
// 许可证的约束，该许可证可以在许可证文件中找到。

package math

// 下面的原始C代码、长注释和常量
// 来自FreeBSD的/usr/src/lib/msun/src/s_log1p.C 
// 并随此通知而来。go代码是原始C.
// 
// ================================================================================================
// 版权所有（C）1993由Sun Microsystems，Inc.保留所有权利。
// 
// 由SunPro开发，SunPro是Sun Microsystems，Inc.的一家企业。
// 使用、复制、修改和分发本
// 软件的权限是免费授予的，前提是保留本通知
// 的权利。第19章
// 1+x=2**k*（1+f），
// ，其中sqrt（2）/2<1+f<sqrt（2）。
// 
// 注。如果k=0，那么f=x是精确的。然而，如果k=0，则f 
// 可能无法精确表示。在这种情况下，需要修正
// 项。设u=1+x四舍五入。设c=（1+x）-u，然后
// log（1+x）-log（u）~c/u。因此，我们继续计算log（u），
// 并加回校正项c/u。
// （注意：当x>2**53时，只需返回log（x））
// 
// 2。log1p（f）的近似值。
// 设s=f/（2+f）；基于log（1+f）=log（1+s）-log（1-s）
// /=2s+2/3s**3+2/5s**5+…，
// /=2s+s*R 
// 我们在[0,0.1716]上使用了一种特殊的Reme算法为了生成近似R的14次多项式，该多项式近似的最大误差为2**-58.45。在
// 换句话说，
// 2 4 6 8 10 12 14 
// R（z）~Lp1*s+Lp2*s+Lp3*s+Lp4*s+Lp5*s+Lp6*s+Lp7*s 
// （程序中列出了Lp1到Lp7的值）
// /和
// /| 214 |-58.45 
// /| Lp1*s+…+Lp7*s-R（z）<=2 
// /| 
// /注意2s=f-s*f=f-hfsq+s*hfsq，其中hfsq=f*f/2。
// 为了保证日志中的错误低于1ulp，我们通过
// 。中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网。最后，log1p（x）=k*ln2+log1p（f）。
// log1p（f）=f-（hfsq-s*（hfsq+R））计算日志
// /=k*ln2_hi+（f-（hfsq-（s*（hfsq+R）+k*ln2_lo）））
// 这里ln2被分成两个浮点数：
// /ln2_hi+ln2_lo，
// 其中n*ln2_hi总是精确表示| n |<2000。
// 
// 特殊情况：
// log1p（x）为NaN，如果x<1（包括-INF），则为信号；
// log1p（+INF）为+INF；log1p（-1）是带信号的-INF；
// log1p（NaN）是没有信号的NaN。
// 
// 精度：
// 根据误差分析，误差始终小于
// 1 ulp（最后一位的单位）。
// 
// 常数：
// 十六进制值是以下
// 常数的预期值。可以使用十进制值，前提是
// 编译器将十进制值转换为二进制值的精度足以产生所示的十六进制值。
// 
// 注意：假设log（）返回准确答案，以下
// 算法可用于在几个ULP内计算log1p（x）：
// 
// u=1+x；
// 如果（u==1.0）返回x；else 
// 返回日志（u）*（x/（u-1.0））；
// 
// 请参阅HP-15C高级功能手册，第193页。

// Log1p返回1的自然对数加上它的参数x。
// 当x接近零时，它比Log（1+x）更精确。
// 
// 特殊情况为：
// Log1p（+Inf）=+Inf 
// Log1p（±0）=±0 
// Log1p（-1）=-Inf 
// Log1p（x<-1）=NaN 
// Log1p（NaN）=NaN 
func Log1p(x float64) float64 {
	if haveArchLog1p {
		return archLog1p(x)
	}
	return log1p(x)
}

func log1p(x float64) float64 {
	const (
		Sqrt2M1     = 4.142135623730950488017e-01  // Sqrt（2）-1=0x3fda827999fcef34 
		Sqrt2HalfM1 = -2.928932188134524755992e-01 // Sqrt（2）/2-1=0xBFD2BEC3333018866 
		Lp5         = 1.818357216161805012e-01     // /3FC7466496CB03DE 
		Lp6         = 1.531383769920937332e-01     // /3FC39A09D078C69F 
		Lp7         = 1.479819860511658591e-01     // /3FC2F112DF3E5244 
	)

	// /特殊情况
	switch {
	case x < -1 || IsNaN(x): // /包括-Inf 
		return NaN()
	case x == -1:
		return Inf(-1)
	case IsInf(x, 1):
		return Inf(1)
	}

	absx := Abs(x)

	var f float64
	var iu uint64
	k := 1
	if absx < Sqrt2M1 { // /| x |</Sqrt（2）-1 
				return x
			}
			return x - x*x*0.5
		}
			// （Sqrt（2）/2-1）<x<（Sqrt（2）-1）
			k = 0
			f = x
			iu = 1
		}
	}
	var c float64
	if k != 0 {
		var u float64
		if absx < Two53 { // 1<<53 
			u = 1.0 + x
			iu = Float64bits(u)
			k = int((iu >> 52) - 1023)
			// 校正项
			if k > 0 {
				c = 1.0 - (u - x)
			} else {
				c = x - (u - 1.0)
			}
			c /= u
		} else {
			u = x
			iu = Float64bits(u)
			k = int((iu >> 52) - 1023)
			c = 0
		}
		iu &= 0x000fffffffffffff
		if iu < 0x0006a09e667f3bcd { // Sqrt的尾数（2）
			u = Float64frombits(iu | 0x3ff0000000000000) // 标准化u 
		} else {
			k++
			u = Float64frombits(iu | 0x3fe0000000000000) // 标准化u/2 
			iu = (0x0010000000000000 - iu) >> 2
		}
	}
	hfsq := 0.5 * f * f
	var s, R, z float64
		if f == 0 {
			if k == 0 {
				return 0
			}
			c += float64(k) * Ln2Lo
			return float64(k)*Ln2Hi + c
		}
		if k == 0 {
			return f - R
		}
		return float64(k)*Ln2Hi - ((R - (float64(k)*Ln2Lo + c)) - f)
	}
	s = f / (2.0 + f)
	z = s * s
	R = z * (Lp1 + z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))))
	if k == 0 {
		return f - (hfsq - s*(hfsq+R))
	}
	return float64(k)*Ln2Hi - ((hfsq - (s*(hfsq+R) + (float64(k)*Ln2Lo + c))) - f)
}
